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We consider a leaky integrate-and-fire neuron with deterministic subthreshold dynamics and a 
firing threshold that evolves as an Ornstein-Uhlenbeck process. The formulation of this minimal 
model is motivated by the experimentally observed widespread variation of neural firing thresholds. 

We show numerically that the mean first passage time can depend non-monotonically on the noise 
amplitude. For sufficiently large values of the correlation time of the stochastic threshold the mean 
first passage time is maximal for non-vanishing noise. We provide an explanation for this effect by 
analytically transforming the original model into a first passage time problem for Brownian motion. 

This transformation also allows for a perturbative calculation of the first passage time histograms. 

In turn this provides quantitative insights into the mechanisms that lead to the non-monotonic 
behaviour of the mean first passage time. The perturbation expansion is in excellent agreement 
with direct numerical simulations. The approach developed here can be applied to any deterministic 
subthreshold dynamics and any Gauss-Markov processes for the firing threshold. This opens up 
the possibility to incorporate biophysically detailed components into the subthreshold dynamics, 
rendering our approach a powerful framework that sits between traditional integrate-and-fire models 
and complex mechanistic descriptions of neural dynamics. 

PACS numbers: 87.19.11, 87.19.1c, 05.40.-a 


I. INTRODUCTION 

In vivo recordings of the membrane voltage in many 
types of neurons display stereotypical upstrokes known 
as spikes [1]. The timing of these spikes has been shown 
to exhibit large variability. While fluctuations are often 
considered as obscuring biological function, the construc¬ 
tive role of stochasticity in neural information processing 
is now well established las]. 

The random nature of spike generation has attracted 
considerable attention in the field of mathematical neu¬ 
roscience. Among the many models that have been 
suggested to describe neural stochasticity, the class of 
integrate-and-fire (IF) models [5] has been used to great 
effect [6]. At the heart of IF models is the notion of 
neural excitability. A voltage spike is elicited when the 
membrane potential reaches a threshold. Traditionally, 
the threshold is considered as either constant or a given 
function of time that can depend on the spike history, 
while the membrane potential is described by stochastic 
differential equations [zHU]. 

In the present study we investigate the timing of spikes 
that are generated when a deterministic subthreshold 
voltage crosses a stochastic firing threshold. A fluctu¬ 
ating threshold reflects experimental findings that the 
membrane voltage at which a spike is elicited varies. 
This has been demonstrated e.g. in cortical neurons 
and hippocampal pyramidal cells na ng. Numerous 
mechanisms have been suggested that could give rise to 
a variable threshold including adaptation, channel noise 
and dynamic modulation of the axon initial segment pTl - 
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EO]. The physiological relevance of a variable threshold 
was recently demonstrated in cortical neurons where syn¬ 
chrony detection was significantly improved [15]. In [21] 
a leaky IF model coupled to a threshold that evolves 
as a Gaussian process successfully described spiking be¬ 
haviour of regularly firing stellate neurons within the ven¬ 
tral cochlear nucleus. A stochastic threshold also cap¬ 
tures inherent uncertainty in both the detection of spik¬ 
ing thresholds and the spike generation mechanism. The 
former arises from non-standardised methods to deter¬ 
mine threshold values from experimental records [22] and 
the fact that spikes are generated at the axon initial seg¬ 
ment, but are often recorded at the soma. The latter 
results from still incomplete knowledge of the molecular 
components that trigger a spike [20] . 

IF models with a fluctuating threshold have been in¬ 
vestigated in the past. However, the approach presented 
here significantly differs from these studies. Often a dis¬ 
tribution of threshold values has been assumed at each 
time point [231426] . As such, successive threshold val¬ 
ues are uncorrelated, and in principle arbitrarily large 
changes of the threshold can occur due to the unbounded 
support of the chosen threshold distributions. In [27] the 
concept of hazard functions is used to determine opti¬ 
mal parameters of the firing probability. The tested haz¬ 
ard functions have unbounded support, and the analysis 
assumes an inhomogeneous Poisson process as the spike 
generation mechanism [28] . A more general approach was 
taken in [29] where a subthreshold stochastic process in¬ 
tersects with a frozen stochastic threshold giving rise to a 
phase transition in the distribution of first passage times. 

In principle it is always possible to transform an IF 
model with a deterministic firing threshold and stochastic 
subthreshold dynamics with additive noise into a model 
with a deterministic description of the membrane voltage 
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and a fluctuating firing threshold (see e.g. [30]). Indeed 
our model is closely related to a leaky IF model driven 
by white noise and subject to an exponentially decaying 
firing threshold m- The time-dependent firing thresh¬ 
old can be transformed into a constant one giving rise to 
an exponential drive in the subthreshold dynamics m- 
Generally such explicitly time-dependent stochastic sub¬ 
threshold dynamics can only be studied numerically. Per¬ 
turbative expressions for the first and second moment of 
the first passage time (FPT) have been obtained in the 
limit of a weak exponential drive [32]. In contrast, we de¬ 
rive analytical expressions of the FPT distribution which 
do not rely on a small parameter. Our findings provide 
quantitative insights into non-trivial effects such as the 
non-monotonic behaviour of the the mean FPT (MFPT) 
as a function of the noise amplitude. 



II. MODEL AND GOVERNING EQUATIONS 

In the subthreshold regime the membrane voltage v{t) 
obeys the ordinary differential equation of a leaky inte¬ 
grator 


I = -«.. +A (I) 

with a and (3 positive constants. A spike occurs when 
v{t) hits the threshold h{t) from below, at which point v 
is discontinuously reset to a constant = i;(0), i.e. 


lim v{t) = Vy . (2) 

We hence define the firing times by 

Tn = mf{t\v{t) > h{t);t > T^-i). (3) 

The fluctuating threshold h is given by 

h = h^sX{t) (4) 

where h > denotes the mean of the threshold and 
5 > 0 measures the coupling strength to a stochastic 
process X{t). We model X{t) as an Ornstein-Uhlenbeck 
process (OUP) whose dynamics reads as [33] 

dX = -jXdt + VDdW. (5) 

Here, 7 is the inverse of the correlation time of the OUP, 
D is a positive constant, and dW is the increment of 
a standard Wiener process W. We choose the initial 
condition X(0) = 0 such that h{0) = h and reset the 
OUP at each threshold crossing to X(0). A fixed value of 
X(0) renders the above model a renewal proeess. With¬ 
out the discontinuous reset of the OUP, consecutive in¬ 
terspike intervals (ISI) are not independent identically 
distributed and hence do not describe a renewal process 
[34]. Fig. 0 illustrates the generation of ISIs by the de¬ 
terministic subthreshold dynamics of v{t) and the OUP 


FIG. 1: (Golor online) Typical sample path of the OUP 
renewal model. After each spike (diamond), the OUP 
and the trajectory are reset to a fixed value 
X(0) = 'r’(O) = 0. Parameter values are 
/3 = 10, h = 9, (a = 7 = 1, D = 2, e = 1. 


in the renewal regime. Note that the full non-st at ionary 
correlation function of the OUP enters our analysis, i.e. 


{Xit)x{t')) = F 




(6) 


for Tn <t^t'< T^+i. It is worth noting that 5 in Eq. Q 
does not need to be small. As both v and X are reset 
at a threshold crossing such that the threshold 

never reaches negative values although the OUP in gen¬ 
eral takes on all values on the real line. Throughout this 
paper we set = 0. 


III. MEAN FIRST PASSAGE TIMES 


We computed the MFPT (T^+i —T^) as a function of 
7 and 5 when [3/a > h. Angular brackets indicate the 
average of the firing times T^+i — for the renewal OUP. 
For notational convenience, we will denote the MFPT 
by T throughout the paper. In this case firing occurs 
even for a deterministic threshold at h. Our results are 
shown in Fig. [^ For small values of 7, the MFPT first 
increases as a function of 5 before decreasing. As we 
increase 7, the range over which T grows monotonically 
as well as the amplitude both decrease. The amplitude is 
defined as the maximal increase of T with respect to the 
noise-free constant firing time. The latter always exists 
for the parameter regime under investigation and can be 
computed from Eqs. 0-0 with 5 = 0 as 


Tdet = — In 
a 


P _ 

/3 — ah 


H 




(7) 
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Here, H denotes the Heaviside step function, i.e. H{x) = 
1 for X > 0 and zero otherwise. For 7 = 0.5 the MFPT 
already decays monotonically. As the OUP can be sim¬ 
ulated exactly [35] and the solution for the subthresh¬ 
old dynamics is known in closed form, the main error 
for computing ISIs results from the determination of the 
threshold crossing. We employed two approaches for this. 
One is based on linear interpolation [36] , while the other 
makes use of a Brownian bridge m Both methods 
give almost identical results (green and red diamonds in 
Fig.§. We used more than 10^ independent realisations 
for the computation of the MPFT at each data point. 

Another approach to compute the MFPT consists of 
solving the partial differential equation (PDF) 

e^Dd^T - , , ar ar ^ 

which describes the MFPT T = T(i;o, ho) in dependence 
on general initial values of the voltage and the threshold, 
i.e. 'r’(O) = vq and h(0) = ho. In the present study we 
are interested in the case ho = h and 'T’o = 0- Equation 
®^can be derived by applying the Feynman-Kac formula 
|38l|39] to the stochastic system given by Eqs. Q-([^. In 
the derivation we used the stochastic differential equation 
for h, which follows from Eqs. and Ito’s formula 

[33] as 


dh = -j{h- h) dt + e\/DdW. (9) 

We solve Eq. on a two-dimensional rectangular do¬ 
main bounded by the line vo = ho, on which T = 0. We 
use no-flux boundary conditions on the remaining three 
sides of the rectangle. The size of the rectangle is chosen 
such that these three sides are sufficiently far away from 
the point of interest {vr, h) eliminating any impact of the 
no-flux boundary conditions on the solution h) [40]. 

Insights into the non-monotonic behaviour of the 
MFPT can be obtained by introducing a new stochas¬ 
tic variable g = e^^X, where the dynamics for X is given 
by Eq. (§. Using Ito’s formula, the dynamics of g obeys 

dg = \/De^^dW, ( 10 ) 


and the time-dependent variance of g follows as 


( 11 ) 


We can interpret Eq. ( 11 ) as the definition of a new time 


s in terms of the original time t. We then define a new 
stochastic process V{s) = g{t{s)) — g{0) in time s where 
t{s) denotes the inverse of s(t), which always exists since 
s{t) is strictly monotonically increasing. The variance of 
V is (U^(s)) — (U(s))^ = s, and together with U(0) = 0 
this renders V a standard Brownian motion. This trans¬ 
formation is discussed extensively in [41] and was first 
used in m- From the threshold crossing condition in 
the new time 5 , = h{t{s)) = h^eX{t{s)), we find 





FIG. 2 : (Color online) MFPT T as a function of e for 
7 = 0.1 (A), 0.3 (B), 0.5 (C) and 1.0 (D). Threshold 
crossings are determined by linear interpolation (green) 
or using a Brownian bridge (red). Black diamonds 
indicate numerical solutions to Eq. The blue dashed 
line denotes the deterministic MFPT given by Eq. 0 . 
Parameter values read as a = 1, /^ = 10, h = 9, D = 2. 


where we have used 

X(s)=e-^*(*)(E(s)+X(0)) , (13) 

which follows from the definitions of g and V. Therefore, 
the crossing of a voltage v(t) through a threshold h(t) 
that is described by an OUP is equivalent to the crossing 
of a new subthreshold voltage v(^ through a standard 
Brownian motion V(s). Figure shows the new sub- 
threhold voltage v(s) for different values of 7 and 5 . We 
observe that for a fixed value of 7 the zero crossing of v 
is independent of 5 . It can be shown that this crossing 
occurs at a value so that corresponds to the determinis¬ 
tic first passage time (Eq. 0 ). Note that the value of so 
grows as we increase 7 . For a fixed s < sq the value of 
v(s) increases as we increase 5 . Therefore a given reali¬ 
sation of the standard Brownian motion V will intersect 
with the subthreshold voltage v earlier for larger 5 than 
for smaller 5 for times 5 < sq . For 7 = 1 almost all cross¬ 
ings of V and V occur before sq, which is reflected in the 
monotonically decreasing behaviour of the MFPT T as a 
function of 5 (Fig. [^. For smaller values of 7 there is a 
significant probability that no threshold crossing occurs 
before sq. Since v(s) is a decreasing function of 5 for fixed 
5 > 5o, a given realisation of V that has not intersected 
with V before so will cross v earlier for smaller values of e 
than for larger values. This realisation of V has a larger 
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FIG. 3: (Color online) Subthreshold voltage v{s) for 
different values of 7 and 5 . Dashed lines correspond to 
7 = 1 and solid lines to 7 = 0.1 for 5 = 0.5 (black), 1 
(green) and 1.5 (red). The inset shows a blow-up for 
small times s. Parameter values are h{0) = h, a = 1, 
^ = 10, h = 9 and D = 2. 


FPT for growing values of 5 . The above arguments hold 
for single realisations of V and hence individual FPT. 
To elucidate the increase of the MFPT for growing 5 we 
need to differentiate between different processes. On the 
one hand more realisations of V could intersect with v 
at times larger than sq. On the other hand the spread 
of FPTs that are larger than 5o grows. To determine 
which factor — or the combination of both — is relevant 
for the current observations we compute numerically the 
FPT distributions. 


IV. DISTRIBUTION OF FIRST PASSAGE 
TIMES 


Our calculation of FPT distributions is based on an 
approach developed by Durbin and Williams [43]. They 
derived an alternating series for the FPT density p of a 
Brownian motion through a curved boundary. Using the 
notation of Sec. [HI] for the transformed dynamics in time 
5 , we have 


p(5) = y](-irift(5) + (-7r,(s), (14) 


i=l 


with 


®(s) = / [ ■■■ [ Ki{s)fi{s)dsi-i---dsi. (15) 

Jo Jo Jo 

Here, we introduced 

■d?y(sj_i) ^sj-O-^SjY 


Ki{s)= n 

j=l 


ds 


Sj-I Sj 


(16) 


where we formally set F(s^) = = 0. fi{s) denotes the 

joint probability density of U(s), U(si),..., U(si_i) at 
the boundary values F(s), F(si),... ,F(5^_i). In Eq. (0 
the term rk{s) represents the error that is made by trun¬ 
cating the infinite series after k terms. Note that Eq. (16) 
differs from the original expression since we study down- 
crossings of a Brownian motion instead of up-crossings 
as in m- Eor practical purposes it is convenient to in¬ 


troduce J^^{s) = We can transform 

results obtained in the new time s to the original time 
t via = J^^{s{t))ds{t)/dt. In Eig. we compare 

with direct numerical simulations for different val¬ 
ues of 7 and 5 . Eor 7 = 0.1 we observe that for increasing 
5 the maximum of the distribution slightly shifts to the 
left, while the tail of the distribution becomes signifi¬ 
cantly longer. This wider spread of long EPTs is the rea¬ 
son for the non-monotonic behaviour of the MEPT as it 
dominates the increase of trajectories with shorter EPTs. 
Eor larger values of 7 the tails of the distributions still 



EIG. 4: (Color online) EPT histograms for 7 = 0.1 (A), 
0.3 (B), 0.5 (C) and 1.0 (D) and e = 0.1 (red), 0.5 
(green) and 1 (blue). The black lines indicates the 
deterministic EPT. The dashed lines correspond to 
J^^(t), while the solid lines result from direct numerical 
simulations. Parameter values are a = 1.0, (3 = 10, 
h = 9, D = 2. 


grow as we increase 5 , but the shift of the distribution to 
shorter times is more pronounced. The latter leads to the 
monotonic decrease of the MEPT for increasing 5 . The 
occurrence of long EPTs for small values of 7 can be at¬ 
tributed to the correlation time of the OUP tqu = 1/7- 
Eor small values of 7, changes in the OUP occur slowly. 
Since for times larger than l/a, the voltage has almost 
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reached its equilibrium /3/(a, the OUP can then spend 
a considerable time in the vicinity of the voltage with¬ 
out crossing it. It is worth noting that captures 

the full FPT distribution extremely well for most of the 
probability mass. only fails at correctly predicting 

the tail of the FPT distribution. The extent to which 
agrees with the true distribution depends on the 
k. As we illustrate in Fig. approximates results 

from direct simulations better as we increase k. In their 
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FIG. 5: (Color online) FPT distributions for 

k = 1 (red), 2 (green) and 3 (blue). The grey line 
corresponds to results from direct numerical 
simulations. Parameter values are 7 = 0.5, e = 0.5, 

D = 2, h = 9, = 10. 

original contribution Durbin and Williams prove conver¬ 
gence of Eq. ( p!4| ) under the assumption that the inter¬ 
cept at the origin of the tangent line of v{s) for any s, 
i.e. v{s)/s — dF(s)/d 5 , does not change sign. Inspection 
of Fig. shows, however, that the intercept is negative 
for small values of s and positive for large values. Our 
results suggest that convergence of the FPT distribution 
also occurs if the above assumption is violated, but a 
rigorous proof is still missing. 

V. CUMULATIVE DENSITY FUNCTIONS 

We next studied the emergence of long tails in the FPT 
distributions in more detail. In Fig. we show the cumu¬ 
lative distribution function (CDF) of the FPTs for differ¬ 
ent values of 7 . We employed the method developed by 
Wang and Potzelberger [44]. The CDF is expressed as 
an expectation value {g) using a piecewise linear approx¬ 
imation of the boundary v. An advantage of this method 
is that each instance of g is known in closed form and 
only the average needs to be computed numerically. We 
observe that for times smaller than 2 the two CDFs agree 
well. This entails that the FPT distributions almost co¬ 



FIG. 6 : (Color online) Cumulative distribution function 
(CDF) of FPTs for 7 = 0.1 (red) and 7 = 1.0 (blue). 
Diamonds represent results using the quasi-analytical 
approach in [44], while solid lines correspond to results 
from direct numerical simulations. Parameter values are 
a = l, e = 1.0, D = 2, h = 9, ^ = 10. 


incide since they follow from the corresponding CDFs by 
differentiation with respect to time. The CDF for 7 = 1 
levels off at significantly earlier times than the CDF for 
7 = 0.1 illustrating the more pronounced tail of the FPT 
distribution for 7 = 0.1. For comparison we also in¬ 
clude results from direct numerical simulations, which lie 
on top of the quasi-analytical results for the CDF. The 
findings in Fig. [^ provide an independent test for the 
non-monotonic behaviour of the MFPT. Moreover, con¬ 
vergence of the CDF as we increase the number of linear 
segments to approximate v is guaranteed. 

The shift of the FPT distributions to smaller times can 
also be understood by computing the probability c that 
a given realisation of V(s) intersects with F(5) at times 
s < 50 - We computed c based on the approach in [44] 
and plot results as a function of 5 for different values of 7 
in Fig.[^ We observe that c is an increasing function of e. 
Hence a growing number of realisations of V have FPTs 
that are smaller than sq. As we increase 7 for fixed 5 we 
find that c grows monotonically demonstrating that more 
realisations of V reach v before sq. The rate of increase of 
c with 5 is more pronounced for larger values of 7 . There¬ 
fore, the relative increase in realisations of V that have 
FPTs smaller than sq is stronger for larger values of 7 
compared with smaller values. Taken together the results 
in Figs. [^ and [^demonstrate the subtle interplay between 
the emergence of long tails in the FPT distribution and 
the increase in the probability mass of FPTs smaller than 
Tdet in the generation of non-monotonic MFPTs. 
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FIG. 7: (Color online) Probability c computed using 
[44] that V(s) has crossed v{s) at least once in the time 
interval (0, sq] as a function of e for 7 = 1 (red), 0.5 
(green), 0.3 (blue) and 0.1 (purple). Parameter values 
are a = 1, D = 2, h = 9, /3 = 10. 


VI. DISCUSSION 

In this paper, we have formulated and investigated an 
IF model for neural firing with a stochastic process con¬ 
trolling the firing threshold, motivated by experimental 
studies showing apparent variability or randomness in the 
spiking threshold. The model is designed to be as simple 
as possible while allowing fluctuations of the threshold 
around a ’preferred’ value h. We therefore use an OUP 
reverting to a mean h for the threshold, and a leaky IF 
model for the sub-threshold dynamics v. The simplicity 
of the model permits the use of techniques from stochas¬ 
tic calculus [38] . 

We have demonstrated that the MFPT may depend 
non-monotonically on the noise amplitude e, finding a 
maximum of the MFPT as the noise strength is in¬ 
creased. The occurrence of a minimum in spiking ac¬ 
tivity at non-vanishing noise strength has been termed 
inverse stoehastie resonanee and has been described in 
stochastic Hodgkin-Huxley models [45]. The emergence 
of inverse stochastic resonance in these models is directly 
linked to the phase-space structure of the correspond¬ 
ing deterministic models. In contrast the deterministic 
LIF model does not predict the complex behaviour of its 
stochastic counterpart. Inverse stochastic resonance has 
been observed in vivo [46] and may even have therapeutic 
applications m- 

Local extrema of the MFPT have been reported be¬ 
fore. In [32| escape from a linear and parabolic potential 
in the presence of a weak exponential driving force is 
considered, while a Kramers problem with a fluctuating 
barrier is investigated in [48] . In both studies extrema of 
the MFPT emerge when a time scale is varied, which is 


the decay rate of the weak exponential forcing in [32| and 
the correlation time of the random barrier in [48|. The 
dependence of the MFPT as a function of the noise ampli¬ 
tude is always monotonic in these studies. This sets our 
results apart from these earlier findings. Moreover, trans¬ 
forming our model to one with a constant threshold and 
an exponential drive in the subthreshold voltage results 
in an amplitude of the time-dependent drive that is not 
small for the parameter values that we have considered. 
This precludes us from applying the perturbative results 
in [SI [32 for computing the first two moments of the 
MFPT. An exponential drive of sufficient amplitude and 
decay time is essential for observing the non-monotonic 
behaviour of the MFPT. The non-mo notonicity results 
from the competition between small passage times, which 
arise when the OUP crosses v during the rising phase of 
the subthreshold voltage, and large passage times, which 
occur when the OUP intersects with the almost station¬ 
ary value of V. 

The dependence of the MFPT on the noise strength 
was found both by direct Monte Carlo simulations and 
by solving a backward Kolmogorov-type PDF, with vir¬ 
tually identical results. The increase of the MFPT with 
noise strength can be explained quantitatively by inspec¬ 
tion of the FPT distributions. As the noise becomes 
stronger, the peak of the FPT distribution moves to¬ 
wards shorter firing times, but the tail of the distribution 
strengthens significantly, so that rare long times between 
firings become more frequent. Further understanding of 
this effect was achieved by transforming the dynamics 
of the system to that of a standard Brownian motion. 
This allows a series method [43] to be used to obtain ap¬ 
proximations to the FPT distribution analytically; this 
series appears to converge rapidly, showing good agree¬ 
ment with numerical results. 

A number of previous studies on LIF models with con¬ 
stant threshold employed different noise intensities of the 
OUP. For example, a noise scaling of was used 

in m, which entails that the variance of the OUP is 
given by D and hence does not depend on the correla¬ 
tion time 1 / 7 . By setting the noise intensity to V 2 D 7 
the authors of [50] ensure that the variance of the sub¬ 
threshold voltage remains almost constant in the limit of 
long correlation times of the OUP, for which the variance 
is given by D 7 . Since the emergence of the maximum of 
the MFPT depends on the interplay of long correlation 
times and sufficient variance of the OUP we tested both 
of the above noise intensities. Our results show that the 
new noise scalings of the OUP have no effect on the non¬ 
monotonicity of the MFPT. The maximum still exists 
but is shifted to larger values of 5 for a given value of 7 . 

An interesting question is how broadly applicable the 
transformation of a stochastic threshold to Brownian mo¬ 
tion is, particularly since an IF model with additive noise 
subthreshold dynamics and constant firing threshold can 
always be transformed into a model where the subthresh¬ 
old voltage evolves deterministically and hits a stochastic 
threshold. If the stochastic threshold is given by a Gauss- 
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Markov process a transformation to Brownian motion al¬ 
ways exists m- A wider class of practically relevant 
stochastic threshold processes are general diffusion pro¬ 
cesses [33] . Here a transformation to Brownian motion is 
only feasible for specific forms of the drift and diffusion 
functions. Wang and Potzelberger m provide a practi¬ 
cal PDE that these functions need to satisfy to derive an 
equivalent Brownian motion. For the most general con¬ 
ditions to recast a diffusion process as Brownian motion 
we refer the reader to [52] . 

Future work will consider two variations of the set-up 
of the problem investigated here. We have considered a 
‘renewal’ model in which the threshold is reset to a con¬ 
stant value after firing, but an alternative procedure is 
to set it to a value selected from the statistical distri¬ 
bution of the OUP. A second possibility is to use a full 


spike train, in which the threshold OUP is not reset after 
each spike, but evolves freely. This will give rise to se¬ 
rial correlations of the interspike intervals, and it will be 
of interest to investigate whether the serial correlations 
display stronger or weaker correlations than the OUP 
threshold process itself. 
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